/*
 * LU.cs
 * 
 * Copyright (c) 2003-2005, dnAnalytics. All rights reserved.
*/

#region Using Directives
using System;
using System.Collections.Generic;
using System.Reflection;
using System.Text;
using dnAnalytics.Math;
using dnAnalytics.Exceptions;
using dnAnalytics.Resources;
#endregion Using Directives

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    /// <summary>
    /// Class for computing the <b>LU</b> factorization of an n by n <see cref="Matrix"/>.
    /// </summary>
    public sealed class LU : Algorithm
    {
        #region Fields
        private ILU decomp;
        private readonly int order;
        #endregion Fields

        #region Constructor
        /// <summary>
        /// Constructs an <b>LU</b> object for the given matrix.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
        /// <exception cref="NotSquareMatrixException">If <paramref name="matrix"/> is not a square matrix.</exception>
        public LU(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix", Strings.NullParameterException);
            }

            if (matrix.Rows != matrix.Columns)
            {
                throw new NotSquareMatrixException(Strings.MustBeSquare);
            }

            decomp = new DenseLU(matrix);

            order = matrix.Rows;
        }
        #endregion Constructor

        #region Public Members

        /// <summary>
        /// Returns a value indicating whether the matrix is singular.
        /// </summary>
        /// <returns><b>true</b> if the matrix is singular; otherwise, <b>false</b>.</returns>
        public bool IsSingular()
        {
            return decomp.IsSingular();
        }

        /// <summary>
        /// Returns the lower factor of the factorization.
        /// </summary>
        /// <returns>The lower factor of the factorization.</returns>
        public Matrix LowerFactor()
        {
            return decomp.LowerFactor();
        }

        /// <summary>
        /// Places the lower factor in the result <see cref="Matrix"/>.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to store the lower factor in.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> doesn't have conforming dimensions.</exception>       
        public void LowerFactor(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }
            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            decomp.LowerFactor(result);
        }

        /// <summary>
        /// Returns the upper factor of the factorization.
        /// </summary>
        /// <returns>The upper factor of the factorization.</returns>
        public Matrix UpperFactor()
        {
            return decomp.UpperFactor();
        }

        /// <summary>
        /// Places the upper factor in the result <see cref="Matrix"/>.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to store the lower factor in.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> doesn't have conforming dimensions.</exception>
        public void UpperFactor(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }
            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            decomp.UpperFactor(result);
        }

        /// <summary>
        /// The determinant of the matrix.
        /// </summary>
        /// <returns>The determinant of the matrix.</returns>
        public double Determinant()
        {
            return decomp.Determinant();
        }


        /// <summary>
        /// Estimates the reciprocal the L1 norm condition number of the matrix. The number is calculated 
        /// as ||A|| * ||AInverse|| where ||.|| is the L1 norm.
        /// </summary>
        /// <returns>The reciprocal of the L1 norm condition number (estimated) of the matrix.</returns>
        public double L1NormConditionNumber()
        {
            return decomp.L1NormConditionNumber();
        }

        /// <summary>
        /// Estimates the reciprocal the infinity norm condition number of the matrix. The number is calculated 
        /// as ||A|| * ||AInverse|| where ||.|| is the infinity norm.
        /// </summary>
        /// <returns>The reciprocal of the infinity norm condition number (estimated) of the matrix.</returns>
        public double InfinityNormConditionNumber()
        {
            return decomp.InfinityNormConditionNumber();
        }

        /// <summary>
        /// Returns the inverse of the matrix.
        /// </summary>
        /// <returns>The inverse of the matrix.</returns>
        /// <exception cref="SingularMatrixException">If the matrix is singular.</exception>
        public Matrix Inverse()
        {
            if (IsSingular())
            {
                throw new SingularMatrixException();
            }
            return decomp.Inverse();
        }

        /// <summary>
        /// Places the inverse of the matrix in to the result matrix.
        /// </summary>
        /// <param name="result">The matrix to place the inverse into.</param>
        /// <exception cref="SingularMatrixException">If the matrix is singular.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> doesn't have conforming dimensions.</exception>
        public void Inverse(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }
            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }
            decomp.Inverse(result);
        }


        /// <summary>
        /// Returns an <b>int</b> array indicating which rows were interchanged during factorization.
        /// Row i was interchanged with row pivots[i].</summary>
        /// <value>An <b>int</b> array indicating which rows were interchanged during factorization.</value>
        public int[] Pivots()
        {
            return decomp.Pivots();
        }

        /// <summary>
        /// Solves a system of linear equations, <b>AX = B</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <returns>The left hand side <see cref="Matrix"/>, <b>X</b>.</returns>
        /// <exception cref="SingularMatrixException">If the matrix is singular.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="rightHandSide"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If the number of rows of the <paramref name="rightHandSide"/> 
        /// matrix differs from the number of rows of the decomposed matrix.</exception>
        public Matrix Solve(Matrix rightHandSide)
        {
            if (rightHandSide == null)
            {
                throw new ArgumentNullException("rightHandSide", Strings.NullParameterException);
            }

            if (rightHandSide.Rows != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            return decomp.Solve(rightHandSide);
        }

        /// <summary>
        /// Solves a system of linear equations, <b>AX = B</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
        /// <exception cref="SingularMatrixException">If the matrix is singular.</exception>
        /// <exception cref="ArgumentNullException">If either <paramref name="rightHandSide"/> or <paramref name="result"/> 
        /// is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If the number of rows of the <paramref name="rightHandSide"/> 
        /// matrix differs from the number of rows of the decomposed matrix.</para>
        /// <para> - or - </para>
        /// <para>If the <paramref name="rightHandSide"/> matrix has different dimensions than
        /// the <paramref name="result"/> matrix.</para>
        /// </exception>
        public void Solve(Matrix rightHandSide, Matrix result)
        {
            if (rightHandSide == null)
            {
                throw new ArgumentNullException("rightHandSide", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (rightHandSide.Rows != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (result.Rows != this.order || result.Columns != result.Columns)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            decomp.Solve(rightHandSide, result);
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side vector, <b>b</b>.</param>
        /// <returns>The left hand side <see cref="Vector"/>, <b>x</b>.</returns>
        /// <exception cref="SingularMatrixException">If the matrix is singular.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="rightHandSide"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If the size of <paramref name="rightHandSide"/> does not equal
        /// the number of rows of the decomposed matrix.</exception>
        public Vector Solve(Vector rightHandSide)
        {
            if (rightHandSide == null)
            {
                throw new ArgumentNullException("rightHandSide", Strings.NullParameterException);
            }

            if (rightHandSide.Count != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            return decomp.Solve(rightHandSide);
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Vector"/>, <b>b</b>.</param>
        /// <param name="result">The left hand side <see cref="Vector"/>, <b>x</b>.</param>
        /// <exception cref="SingularMatrixException">If the matrix is singular.</exception>
        /// <exception cref="ArgumentNullException">If either <paramref name="rightHandSide"/> or <paramref name="result"/>
        /// is <b>null</b>.</exception>
        /// <exception cref="NotConformableException"><para>If the size of <paramref name="rightHandSide"/> does not equal
        /// the number of rows of the decomposed matrix.</para>
        /// <para> - or - </para>
        /// <para>If the size of <paramref name="result"/> does not equal
        /// the number of columns of the decomposed matrix.</para>
        /// </exception>
        public void Solve(Vector rightHandSide, Vector result)
        {
            if (rightHandSide == null)
            {
                throw new ArgumentNullException("rightHandSide", Strings.NullParameterException);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (rightHandSide.Count != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (result.Count != this.order)
            {
                throw new NotConformableException(Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            decomp.Solve(rightHandSide, result);
        }

        #endregion Public Members

        #region Protected Members

        /// <summary>
        /// Delegates the actual computation to the correct LU subtype.
        /// </summary>
        protected override void InternalCompute()
        {
            decomp.Compute();
        }
        #endregion Protected Members

    }
}
